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A b stract 

A four-node, quadrilateral smoothing element is 
developed based upon a penalized-discrete-least-squares 
variational formulation. The smoothing methodology 
recovers C ^continuous stresses, thus enabling effective 
a posteriori error estimation and automatic adaptive 
mesh refinement. The element formulation is originated 
with a five-node macro-element configuration consisting 
of four triangular anisoparametric smoothing elements 
in a cross-diagonal pattern. This element pattern enables 
a convenient closed-form solution for the degrees of 
freedom of the interior node, resulting from enforcing 
explicitly a set of natural edge- wise penalty constraints. 
The degree-of-freedom reduction scheme leads to a very 
efficient formulation of a four-node quadrilateral 
smoothing element without any compromise in robust- 
ness and accuracy of the smoothing analysis. The appli- 
cation examples include stress recovery and error estima- 
tion in adaptive mesh refinement solutions for an 
elasticity problem and an aerospace structural 
component. 

In trodu ction 

The displacement-based finite element method has 
become a standard tool for structural analysis. However, 
due to its inherent interelement stress discontinuity, it 
generally requires a special post-processing or ‘recovery’ 
procedure to enable improved stress calculations. The 
improved stress predictions can further be utilized in a 
posteriori error estimation that in turn enables 
automatic mesh adaptivity. The recovery procedures 
developed to date can be classified as (a) local (i.e., 
element-level), (b) patch-based, or (c) global. The 
simple averaging of nodal stresses from adjacent 
elements is an example of the simplest patch-based 
stress recovery. Early attempts to employ global least- 
squares procedures to obtain continuous stress fields 
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were not widely adopted. 1,2 Instead, local element-level 
least-squares approaches appeared to be more attractive 
because of their intrinsic computational efficiency. 2,3 
However, the local schemes do not recover a continuous 
stress field, and a subsequent nodal averaging is 
generally used to achieve the stress continuity along 
element interfaces. With the recognition that certain 
interior stresses exhibit superconvergent properties, 4 
there evolved a class of local procedures involving 
‘projection’ of these stresses to the element boundary, 
although other local procedures were also proposed. 5 

Recent developments in stress recovery procedures 
have largely been motivated by a posteriori error 
estimation, 6 especially for use in automatic adaptive 
mesh refinement. In this context, polynomial patch 
recovery procedures have been explored. Zienkiewicz and 
Zhu 7,8 developed the Superconvergent Patch Recovery 
(SPR) procedure, and various modifications of their 
approach have been proposed subsequently. 9 ' 11 These 
procedures generally attempt to recover from the 
superconvergent stress points a stress field with 
superconvergent properties. The latest formulation by 
Boroomand and Zienkiewicz 12,13 explored yet another 
polynomial patch-based procedure which uses nodal 
equilibrium rather than superconvergent stresses. In 
general, patch procedures are used to recover the stresses 
at nodes. If a node belongs to multiple patches, an 
averaging procedure must be used. The continuous 
stress field is then defined in an ad hoc manner by the 
nodal stresses and the finite element displacement 
interpolation functions. 

A major departure from other procedures is the finite 
element-based recovery methodology of Tessler and co- 
workers. 14 ' 21 The methodology, which is based on the 
minimization of a Penalized-Discrete-Least-Squares 
(PDLS) error functional, is highly effective and is 
designed to provide improved stress predictions with a 
higher degree of smoothness and to obtain robust a 
posteriori error estimates. 

The error functional involves a discrete least-squares 
term in which discrete finite element stresses are 
compared with continuous recovered stresses, a penalty 
constraint term that enforces C 1 -continuity of the recov- 
ered stresses, and a curvature-control term that ensures 
stability and robustness of the method. The post- 
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processing is performed via the Smoothing Element 
Analysis (SEA) which itself is a finite element method. 
As such, SEA has its own mesh and element-based 
interpolation functions which minimize the PDLS error 
functional. The methodology is capable of recovering 
superconvergent stresses of higher accuracy and continu- 
ity than the original finite element stresses. The 
approach is unique in that the recovered stress field is 
essentially C 1 -continuous. In addition to achieving a 
higher degree of continuity, SEA generally produces 
more accurate stresses than the SPR method. The 
smoothed stresses can be used for design since they 
represent more accurate stresses than the ones obtained 
directly from a finite element analysis. They can also be 
employed in a posteriori error estimators. Importantly, 
SEA fits perfectly in any general purpose finite element 
code since it has the same architecture as the standard 
finite element method. 22 ' 24 

The smoothing element developed by Tessler and co- 
workers is a three-node triangle, incorporating a 
quadratic interpolation for the stress, and linear 
interpolations for the two normal stress derivatives. 20,21 
Since the PDLS error functional possesses a penalty- 
constraint term, the use of anisoparametric 
interpolations is key to enforcing pointwise C 1 stress 
continuity without developing the detrimental locking 
effect. The method is especially robust and effective 
when the smoothing element discretization is 
represented by quadrilaterals consisting of four triangular 
smoothing elements in a cross-diagonal pattern. The 
resulting macro-element consists of five nodes, each 
having three degrees-of-freedom (dof s). 

In this paper, a computationally efficient, four-node 
quadrilateral smoothing element is developed. The 
element is derived from a five-node macro-element 
consisting of four triangular smoothing elements in a 
cross-diagonal pattern. This element pattern enables an 
exact closed-form solution for the dof s of the interior 
node, resulting from enforcing explicitly a set of natural 
edge- wise penalty constraints. The degree-of-freedom 
reduction scheme gives rise to a four-node quadrilateral 
smoothing element which has obvious computational 
advantages over the five-node macro-element. In 
particular, its mesh results in nearly one-half the dof s 
than would a corresponding cross-diagonal mesh 
constructed with triangular elements. Because an exact 
closed-form solution is used for the reduction, no 
compromise is made in terms of robustness and 
accuracy by reducing the dof s. 

The conceptual framework and variational basis of 
SEA/PDLS are first discussed. Then, smoothing equa- 
tions for a triangular element are derived on the basis of 
consistent anisoparametric interpolations for the stress 
and two normal stress derivatives. Eight edge- wise 


constraint equations are formulated for a five-node 
macro-element consisting of four triangles in a cross- 
diagonal pattern. Three of these equations are solved in 
closed form for the dof s corresponding to the interior 
node. The resulting transformation relations are used to 
formulate the four-node smoothing element in a 
computationally efficient manner. Numerical examples 
include stress recovery and error estimation in adaptive 
mesh refinement solutions for an elasticity problem and 
an aerospace structural component. 

SEA Con cep tu a 1 F ram ew o rk 

The basic conceptual framework of the present stress 
recovery and error estimation approach is illustrated in 
Fig. 1. Based upon a Finite Element Analysis (FEA), 
which can be quite general and involve either linear or 
nonlinear solution procedures, a subsequent Smoothing 
Element Analysis (SEA) is carried out which discretizes 
the same physical domain of the problem. The SEA 
employs the discrete FEA stresses c h as the prescribed 
input data, using either Gauss integration point stresses 
or superconvergent stresses. These stress locations are 
preferred because of their superior accuracy compared to 
other stresses. Alternatively, quantities other than stress 
components can also be smoothed in the same manner, 
for example, discrete displacement data, 14 strain 
measures, 19 stress invariants, element strain energies, or 
any other discrete data of interest. 

The SEA itself is a finite element method and 
requires its own mesh which in general can differ from 
the finite element analysis mesh. Thus, the problem 
domain Q is discretized with n el smoothing finite 

elements such that £2 = u^j £2 e , where £2 e is the 

domain of a smoothing element. Within each smooth- 
ing element, a penalized-discrete-least-squares (PDLS) 
error functional is formulated and is based upon a 
smooth stress G s and two stress-gradient variables. 
These field variables are interpolated with compatible 
shape functions that ensure a practically C 1 -continuous 
stress field and C° -continuous stress gradients. The 
total PDLS error functional contributed by all smooth- 
ing elements in the discretization is minimized with 
respect to all dof s of the smoothing field variables. 
The resulting smoothing equations, which have a stan- 
dard linear algebraic form, produce an improved and 
pointwise smooth stress field. The recovered smooth 
stress field can either be used directly for design calcula- 
tions or they can be employed for a posteriori error 
estimation in an automatic adaptive refinement process, 
which is perhaps of greater significance. 
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Error Estimation for Adaptive Refinement 

- Compute error norms n s h n 

- Generate adaptive mesh F ~~ a II 


Fig. 1. Smoothing element analysis for stress recovery and error estimation. 


In developing a practical and theoretically appropriate 
framework for the application of SEA to built-up 
structural models, such as those found in aerospace and 
automotive structures, it is important to realize that 
stresses belonging to different structural components 
may, in fact, be physically discontinuous where such 
components intersect. Moreover, in most cases, stress 
components in one shell or beam segment would lack 
the same definition in another segment. For these rea- 
sons, it may not be meaningful to employ a scheme 
which would require continuity of stresses along struc- 
tural interfaces. Also, the usual averaging of stresses 
along structural interfaces would often result in errone- 
ous stress results at the given interface and the surround- 
ing region. 

The application of SEA to built-up structures is 
undertaken herein by associating individual structural 
segments with an independent SEA mesh. This 
approach should be applicable to one-dimensional 
(beams and trusses), two-dimensional (plates and 
shells), and three-dimensional (solid element models) 
structural finite element approximations. Thus the 
SEA discretization would consist of distinct domain- 
based meshes, each totally independent of the adjacent 
SEA domain. The SEA domain-based framework will 
be further illustrated on a built-up aerospace structure 
discussed in the Numerical Results section. 

E rror Functional 

The mathematical foundation of SEA is described for 
a stress field defined on a flat two-dimensional domain, 


£1 = {x E 9t 2 }, for which x = {x,y} denotes the 
position vector in Cartesian coordinates. As will be 
shown subsequently, this two-dimensional framework 
can be readily extended for application to typical built- 
up aerospace shell structures. 

For a single smoothing element within a SEA 
mesh, the PDLS error functional may be expressed for a 
given stress component G as 

* e =^I [tfq-° S ( X q)f + 

q=l 

n e 

a Ife [(< x -e s x ) 2 +«yf y -e s y ) 2 ]dn + 

e=l 

n e ^ 

PL^ e f e t( 0 lx ) 2 +( 0 y >y ) 2 + 

e=l 

i(e!, y + e s y;X ) 2 ]dQ 

where, for simplicity, the usual tensorial subscripts for 
the stress are omitted; a comma denotes partial 
differentiation, OC and (3 are dimensionless parameters, 
and N denotes the total number of sampled discrete 
stresses. For smoothing purposes, G is treated as a 
scalar quantity. As will be demonstrated shortly, the 
variables 0- (i=x,y) represent the first-order partial 

derivatives of G s with respect to the spatial coordinates. 
The finite element stress field is evaluated (or sampled) 
in each element at x q (q = 1, 2, ..., n e ), to obtain the 
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set of discrete stresses {G q }, i.e., Gq=G h (x q ), and 

where n e is the number of the sampled stresses within 
the element. Because the highest partial derivative in (1) 
is of order one, the field variables need only be 
approximated with C° -continuous shape functions. 
The functional is minimized with respect to all nodal 
dof s and yields a system of linear algebraic equations 
which must be solved for the dof s. 

The first term in (1) is a normalized discrete least- 
squares functional in which the squared ‘error’ between 
the smoothed stress field and the sampled stress data is 
computed for all sampled stresses. The second term in 
(1) represents a penalty functional which, for 
sufficiently large a, enforces the derivatives of the 
smoothed stress field G^ to approach the corresponding 

0- variables pointwise, i.e., 

G^i - >6i (i = x, y) in £2 e . (2) 

The greater the value of a, the closer the correlation 
between G^ and 0-, where C 1 continuity of G s is 

achieved as a — » °o . For the practical application of the 
method, a needs to be sufficiently large in order to 
enforce conditions (2); however, it should not be 
excessively large to cause ill-conditioning of the 
smooth solution. The range a=10 2 — 10 6 has been 
demonstrated in previous studies to be quite adequate for 
these purposes. Because 0- are interpolated with 
continuous functions, the smoothed stress field, for all 
practical purposes, can achieve C 1 continuity 
throughout the Q domain. 

The third term in (1) involves squared first deriva- 
tives of the 0- variables which, in accordance with (2), 
represent the curvatures (i.e., second partial derivatives) 
of the smoothed stress field in the limiting sense, (i.e., 
Gjj — > Oi,j) F° r the shape functions subsequently 

considered for the smoothing element, the 0-j terms 

represent the curvatures of the smoothed stress exactly. 
The third term imposes a constraint condition on the 
stress-curvature field, the severity of which is governed 
by the value of the parameter (3. When the sampled 
stress data is perceived to be reasonably accurate, the (3 
parameter needs to be very small, particularly in relation 
to the penalty parameter a. For the sampled stress data 
exhibiting substantial error, larger values of (3 need to 
be used to smooth or filter the data. A more general 
form of (1) which enables the use of appropriate 
weighting functions for the first and second terms can 
be found in Ref. 20. 


A notable observation regarding the functional form 
of the penalty-constraint and curvature-control terms in 

(1) is that they can be seen as mathematical analogs of 
the transverse shear and bending energies in Reissner- 
Mindlin plate theory. In this connection, successful 
element technologies developed for plate elements can 
equally be applied in this formulation. 25-27 

Sm oothinq E lem ent Interpolations 

The description of a convenient and effective 
interpolation strategy for the field variables to be 
approximated in (1) is addressed. Although the 
functional in (1) admits C° -continuous shape functions 
for the field variables G s , 0| , and 0 y , the constraints in 

(2) severely limit the suitable choice of the shape 

functions. The reasons for these restrictions lie in the 
character of these constraint equations. Unless appro- 
priate interpolations are selected for the field variables, 
the resulting solution will suffer from ‘locking’. In the 
context of the smoothing analysis, ‘locking’ would be 
manifested by the smoothed solution that grossly 
underestimates the stress field and is thus rendered to be 
useless. The locking effect and the computational means 
of avoiding it have been studied extensively, especially 
in the context of shear-deformable bending finite 
elements. 25-27 Although the earlier approaches favored 
reduced integration of the penalty term (a counterpart of 
the second term in (1)), the improved understanding of 
this phenomenon in recent years seems to favor the use 
of special interpolation schemes that ensure consistency 
both in the constraint and variational sense. Such 
interpolation approaches, interchangeably called inter- 
dependent and anisoparametric , were originally 

introduced by Tessler and co-workers 26-28 in their shear- 
deformable beam and plate formulations. 

The three-node triangular smoothing element 
employs the anisoparametric interpolations of the 
lowest order, that is G s is interpolated with a complete 
quadratic polynomial and the 0 X and 0 y employ linear 

interpolations. In reference to the notation in Fig. 2, 
these smoothing element interpolations may be 
expressed in matrix form as : 


g s = Ps + Q x s x +Q y s y = Nd e , 
0| = PS X , 0 y = PS y 


(3) 


where s = [c v g 2 , g 3 ] t , s x = [0 xl , 0 x2 , 0 x3 ] T , and 
s y = [0 yl , 0 y2 , 0 y3 ] T are the vectors of nodal dof’s, and 
P, Q x , and Q y are the row- vectors of linear and 
quadratic shape functions, respectively. The inter- 
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polation functions given in terms of triangular area- 
parametric coordinates, z = [£j, £ 2 , £ 3 ], can be written 
as 


Pi=?i. Qxi =^(a k CiCj -aj^i^k), 

Qyi=^(bjCi^k-b k ^iCj) (4) 


in which 


Ci = 


Cj + bj x + aj y 


2A 


c i = x jy k -x kyj 


ai = x k - Xj, bj =yj -y k , 
A = ^( a 3 b 2- a 2b 3 ), 


with A representing the area of a triangular element, and 
x k , y k the nodal Cartesian coordinates; in the foregoing 
equations, the subscripts are given by the cyclic 
permutation of i = 1, 2, 3, j = 2, 3, 1, and k = 3, 1, 
2. Note that these interpolations are consistent with a 
three-node element which has only three dof s per node. 



with respect to the nodal dof s (i.e., 8d> e =0), giving 
rise to 


K e d e = F e (5) 

where d e is a vector containing the element nodal dof s, 
K e is a symmetric, element smoothing matrix, and F e 
is a consistent right-hand-side vector. The element 
matrices are given as 


and 


K e = Kg + + Kp 


= ^I N q N q + «LeBXd^ + 

^ q=l 

(3n e { BpDBp d£l , 


F e = — V ajj nJ, d e = 

N q tl q q 


L s yj 


Px Qxx-P Qy,x 
Py Qx,y Qy,y- P .’ 


0 

Px 

0 ’ 


1 

0 

o" 

0 

0 

p y 

, D = 

0 

1 

0 

0 

Py 

Px 


0 

0 

1 

2 _ 


(6) 


Fig. 2. Notation for three-node anisoparametric 
smoothing element. 

The following attributes delineate the most notable 
benefits of the three-node anisoparametric triangular 
element for utilization in SEA: (a) The element 
interpolations (3) and (4) ensure that the gradient of the 
smoothed stress, Gj , is the same degree polynomial as 

that representing 0-, i.e., they are both linearly 
distributed across Q e ; this interpolational consistency is 
the key to the proper element-level resolution of 
constraints (2), (b) The triangular element permits a 
one-to-one linear mapping between the global and 
element (area-parametric) coordinates, thus allowing a 
straightforward mapping of the sampled stress data onto 
the smoothing element, and (c) The quadratic 
distribution of G s is relatively low order and thus 
ensures sufficient robustness of the method both in 
terms of interpolation and extrapolation. 

Sm oothincr E quations 

The element smoothing equations are obtained by 
introducing (3) into (1) and taking the first variation 


where the element matrix K e is comprised of the error 
matrix, Kg, the penalty-constraint matrix, K^, and the 
curvature-control matrix, Kp; and N q =N(x q ). Equation 

(6) is integrated using exact integration formulas. Note 
that integration of Kp is trivial since its integrand is 

constant. 

The global smoothing equations, Kd=F, are 
obtained using the usual finite element assembly 
operation and are solved for each stress component. 
Conveniently, K e depends only on the element shape 
functions that are evaluated at the sampled stress 
locations and does not involve the stress values. This 
means that the smoothing of an individual stress 
component corresponds to a different right-hand side 
only, stored in the F vector. Thus, regardless of the 
number of smoothed quantities, the global K matrix 
need only be assembled and factored once. Employing a 
very small (3 parameter (say, 10' 5 * ) assures that K is 
nonsingular even when the minimal number of 
sampling stresses, N=3, is available. 17 Although this 
aspect permits a great range of possibilities for 
constructing a SEA mesh, it is highly desirable for 
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automation purposes that the SEA mesh be identical to 
the FEA mesh. 

F ou rTJ ode Q u adrilateral via E dqe 
C on strain ts 

The three-node triangular interpolations described in 
the previous section are particularly effective and robust 
when the smoothing element discretization is locally 
represented by quadrilaterals consisting of four triangular 
smoothing elements in a cross-diagonal pattern. 20,21 The 
resulting macro-element has five nodes, and a total of 
15 dof’ s. Such a macro-element, in which large penalty 
numbers ensure virtually exact satisfaction of the 
constraints (2), can be used without causing any 
locking-type deterioration, or even degradation in 
accuracy. 

Once (3) are introduced into (2), a straightforward 
manipulation of the resulting equations produces three 
edge- wise constraints per element. It then follows that 
the five-node macro-element (refer to Fig. 3) would 
possess a total of eight such edge constraint equations - 
four interior and four exterior. Replacing the limiting 
condition with that of equality in (2), which for large 
values of a should be valid for all practical purposes, 
the edge constraint equations can be written in the form: 

Interior E dqe C onstraints 


a i ~ a 5 ~ a i5 (®xi +9 x 5) + bi5(9 y i +9y5)’ (?) 

i= 1,2, 3,4 


and 


b = g(yi + y2+y3+y4- 4 y5)- 


Furthermore, a straightforward manipulation of the 
equations in (7) and (8) yields the matrix equation for 
the 0 x5 and 0 y5 dof s, i.e., 




— a 25 

0 

T 

0x1 


a 13 

_ a 35 


0 x2 


a 25 

a 24 


0 x3 

b 24 -b 13 

0 

a 35 


®x4 

— a 24 a 13 

— b 25 

0 


0 yi 


b 13 

— b 35 


0 y2 


b 25 

b 24 


0 v3 


0 

b 35 . 


CD 


where A Q = 2(a 13 b 24 -a 24 b 13 ) is the area of the 
quadrilateral macro-element. Equations (9) and (10) 
provide the explicit relations between the dof s at the 
cross-diagonal node (node 5) and the vertex nodes of the 



Fig. 3. Reduction of five-node macro-element to four- 
node quadrilateral via edge constraints. 


where a i5 =-(x ; -x 5 ), b i5 = -(y ; -y 5 ), 

E x terio r Edge Con stra in ts 

- Gj = ay (0 xi + 0 x j ) + by (0 yi + 0 y j ) , (8) 

i = 1, 2, 3, 4; j = 2, 3, 4, 1 

where ay = ^(x { -xj), by = -^-(yi -yj) 

Adding the four equations in (7) results in the solution 
for g 5 , 

a 5=7 T k - a i5 ©Xi - b i5 e yi ] - a e x5 “ b 0 y5 ( 9 > 

^i=l,4 


where 


a = — (x 1 + x 2 +x 3 +x 4 -4 x 5 ) 


There are two equivalent approaches of employing 
the reduction equations (9) and (10) to construct the 
four-node quadrilateral: (a) The element matrices for the 
four triangles comprising the five-node macro-element 
are pre-imposed and multiplied by the appropriate 
transformation matrices, Tj (j=l,4), resulting from (9) 
and (10). This process is variationally consistent and is 
commonly employed whenever transformations of dof s 
are performed, (b) Alternatively, (9) and (10) may be 
used directly to modify the interpolation functions for 
each of the four triangles in the five-node macro- 
element. Using the foremost method, the quadrilateral 
element equations are as follows: 

Kq dQ = Fq (11) 

where 

K Q= LTfKfTj, F^IT/Ff. 

j=l,4 j=l>4 

and 

d e Q=k’ e *i’ 0 yi} T < i=1 ’ 4 )’ 
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and where the summation is over the four triangular 
elements. 

It should be noted that static condensation of the 
interior dof’ s would result in the same solution as the 
quadrilateral smoothing solution due to (11), however 
the systems of equations are in different bases and are 
therefore not identical. In addition, static condensation 
is not as computationally efficient as the explicit 
reduction equations (9) and (10). 

Num erica 1 E x am p les 

Two adaptive mesh refinement solutions are carried 
out to demonstrate the robustness and accuracy of the 
smoothing method and the computational efficiency of 
the four-node smoothing element. First, a linearly 
elastic plate with a small central hole under 
compression is analyzed to verify the equivalence of the 
five-node, cross-diagonal pattern macro-element and the 
four-node quadrilateral, and to demonstrate the 
computational efficiency of the four-node element. The 
second numerical solution deals with a built-up 
aerospace structural component under compression and 
bending, where shell elements are located in different 
planes. In this structure, stresses at shell junctions are 
in general discontinuous, and for this reason, SEA is 
performed on different domains independently. 

All computations are obtained with NASA's 
COMET-AR (Computational MEchanics Testbed - 
Adaptive Refinement) FEA code on an IBM RS6000 
workstation using nine-node fully integrated Lagrange 
shell elements. 22 These shell elements are somewhat 
stiff in bending, however, they are relatively insensitive 
to mesh distortion - an aspect that is paramount in 
adaptive mesh refinement. No special consideration is 
taken for modeling shell junctions. The smoothing 
mesh generation is fully automatic and uses the same 
quadrilateral mesh as that of the finite element models. 
The FEA stresses are sampled at the 3x3 Gauss points. 
In the smoothing mesh, these points are located with an 
efficient search algorithm. 24 

E rror E stim ation and Adaptive 
R efinem en t 

To demonstrate robustness of the smoothing element 
due to element distortion, aspect-ratio tolerances are not 
set in the adaptive refinement. Instead, the maximum 
number of refinement meshes is specified. 

Adaptive mesh refinement is performed using 
transition-based refinement strategies with a specified 
refinement tolerance. Any finite element exhibiting an 
average element error, R ave , greater than the specified 


refinement tolerance is subdivided into nine elements. 
The average element error is computed as 


Rave= — f (12) 



in which E e is the element error energy norm given as 


E e 


- | (S s - S h ) T C _1 (S s - S h ) d£2 

2 ^e 


(13) 


The vectors S h and S s contain, respectively, eight FEA 
and SEA shell stress components, N elts is the number of 
finite elements, and C is an 8x8 constitutive matrix. In 
(12), U ref is the reference strain energy defined as the 
finite element strain energy corrected by the global 
error 23 


where 


u ref = U h + E 2 


E = 


T E e /U S 


N-i 


(14) 


with U h and U s denoting the total strain energies 
corresponding to the FEA and SEA stress fields, 
respectively. 


P late W ith H ole 

A two-dimensional elasticity problem is analyzed to 
demonstrate the robustness and efficiency of SEA in the 
context of automatic adaptive refinement. A rectangular 
aluminum plate (thickness = 0.1 in., E = 10 Msi, v = 
0.3) is subjected to a tensile uniform displacement of 
0.1 in. along an edge parallel to the y-axis that is 
constrained to move only in the x direction (refer to 
Fig. 4). The opposite edge is clamped and the edges 
parallel to the x-axis are free. The problem is 
particularly challenging because the hole is very small, 
d/w=0.05, where d is the hole diameter. This aspect 
tends to produce distorted elements around the hole 
where a stress concentration takes place. Stress 
singularities at the plate comers, where the displacement 
boundary conditions are prescribed, provide another 
challenge for the adaptive refinement. 21 

The initial mesh (Mesh 0) and three consecutive, 
automatically generated meshes are depicted in Fig. 4. 
At each refinement step, finite elements of greater 
distortion and smaller size are generated. This is 
especially true near the hole and at the plate comers - 
the regions of high stress gradients. By applying SEA 
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Fig. 4. Plate with hole. Initial (Mesh 0) and three refinement meshes. 
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Fig. 5. Percent difference in N x at point A due to five- 
and four-node quadrilateral SEA solutions. 
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Fig. 6. Comparison of CPU time for SEA with four- 
and five-node quadrilaterals. 

for each mesh, robust stress smoothing and error 
estimation is achieved. For Mesh 3, possessing 4,476 
finite element dof’s, a global error norm of less than 1% 
is obtained. 


In Fig. 5, the equivalence between the five- and four- 
node quadrilateral smoothing solutions for Mesh 0 is 
demonstrated. At point A located at the top of the hole 
perimeter, where the axial stress resultant N x is 
maximum, the solutions approach the same value as the 
penalty parameter a becomes larger. Apparently, the 
four-node solution does not suffer from any loss of 
accuracy as a result of the explicit dof reduction inherent 
in the element formulation. Even at relatively small 
values of a (a=10), the two SEA solutions are within 
0.0003%. This excellent degree of correlation is typical 
for the highly stressed regions in the plate. 

Having fewer dof s and achieving equivalent results, 
the four-node smoothing element is significantly more 
efficient than the original five-node macro-element. 
This aspect is clearly demonstrated in Fig. 6 which 
compares CPU time for the two SEA modeling 
approaches. Evidently, the larger the mesh, the greater 
the computational savings attained. Thus, for Mesh 3, 
the SEA computation is reduced by a factor of 6.6 with 
the use of the four-node smoothing quadrilateral 
element. 

H at-S tif fen ed Panel 

The hat-stiffened panel is an example of a built-up 
aerospace structural component. The panel dimensions 
and boundary conditions in Fig. 7 are shown on the 
initial (Mesh 0) discretization. The right edge of the 
panel is subjected to a compressive axial (x-direction) 
displacement of 0.1 in., including the stiffener edge. 
The opposite edges of the panel and stiffener are fixed. 
Along the side edges, symmetric boundary conditions 
are prescribed. Whereas FEA is carried out with a nine- 
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node Lagrange shell element, SEA is performed using 
the four-node quadrilateral with the element parameters 
set at a = 10 5 and (3 = 0. 

From the structural response standpoint, the panel 
has a number of challenging aspects: (a) the stiffness 
eccentricity due to the hat stiffener located on only one 
side of the panel causes local bending, (b) stress 
concentrations occur along the hole perimeter, and (c) 
singular stresses develop at the ends of the stiffener- 
panel junctions, where stresses are transferred according 
to the shear-lag mechanism. It must be noted that the 
initial mesh is very coarse and does not address any of 
the aforementioned issues. 

Knowing that stress fields at the stiffener-panel shell 
junctions are generally discontinuous, the application of 
SEA is carried out on separate, independent domains to 
avoid unphysical smoothing across these junctions. In 
Fig. 8, the panel modeled with Mesh 0 is broken-up 
into nine independent smoothing domains. Within each 


domain, stress resultants are defined to correspond to a 
convenient coordinate system, and these stress compo- 
nents are smoothed independently. At each refinement 
step, error estimation is undertaken based on the domain 
smoothed stresses. The distributions of the average 
element error (refer to (12)) corresponding to the initial 
and two adaptive refinement discretizations are shown in 
Fig. 9. 

An examination of specific stress distributions 
computed via FEA and subsequent SEA for error 
estimation, further validates the effectiveness and 
robustness of the SEA stress recovery. The typical 
stress results are depicted in Figs. 10 and 11, which also 
display the number of dof’ s and CPU time for the entire 
FEA and SEA computations, including model 
definition. The SEA CPU time includes the smoothing 
of all eight shell stress components (i.e., three in-plane 
stress resultants, three bending moments, and two 
transverse shear stress resultants) for each domain. 
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Fig. 9. Hat-stiffened panel mesh. Refinement steps and associated average 

element error. 



Fig. 10. Hat-stiffened panel . M^sbe-ndiitigomomf nt . 


Figure 10 shows the M x bending moment 
distribution, where the FEA results for Mesh 1 are 
averaged at all common element boundaries, enforcing 
stress continuity at the shell junctions. This type of 
post-processing is standard in most general-purpose 
finite element codes. For reference purposes, the results 
for a highly refined mesh (Mesh 2) are used. 


Since the domain-based SEA does not enforce stress 
continuity across domain junctions, the resulting 
smoothed M x distribution is considerably more 
accurate than that of the FEA, strongly resembling the 
reference solution. Observe that SEA-based M x 
distribution is nearly continuous across domains where 
it is expected, without being enforced (refer to the local 
details in Fig. 10). Also note that SEA-based M x 
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Fig. 11. Hat-stiffened panel. Distribution of N xy stress resultant. 



Fig. 12. FEA and SEA (Mesh 1) accuracy comparisons 
for N xy at point A. 

‘uncovers’ the physical stress concentrations that are not 
particularly evident from the FEA results. 

In Fig. 11, the in-plane shear resultant N xy 
distributions for Mesh 1 and Mesh 2 (reference) are 
presented. These results demonstrate once again the 
importance of the independent domain-based smoothing 
for stress post-processing. For example, the N xy stress 
resultant defined in the stiffener which intersects the 
base panel has a different definition than that in the base 
panel. Yet, unless special care is taken, these results 
are commonly averaged at the shell junctions, producing 
erroneous results at the intersections and surrounding 
regions. With the domain-based SEA, however, stress 
continuity is not enforced across the distinct domains 
and physically meaningful stresses are recovered. 
Observe the similarities of the SEA and reference 
distributions. 


An example of quantitative stress improvement 
associated with the SEA recovery is shown in Fig. 12, 
where, corresponding to Mesh 1, N xy is examined in a 
region of high- stress gradient. It is seen that at point A, 
which is located close to the edge of the left stiffener, 
there is a 38% reduction in the error for N xy due to 
smoothing. 



Fig. 13. FEA global percent error in energy norm vs. 
dof’ s. 

It is evident that SEA-based stress recovery is effec- 
tive in identifying stress concentrations in relatively 
coarse FEA models, which allows for more rapid global 
convergence in adaptive mesh refinement. The conver- 
gence of the adaptive refinement process is illustrated in 
Fig. 13. In the figure, the global FEA percent error in 
the energy norm (computed according to (14) by adding 
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all element contributions) is plotted versus the number 
of finite element dof’ s. 

C on clu d in q R em arks 

A four-node quadrilateral smoothing element has 
been developed as a result of an explicit dof reduction 
solution. The original five-node macro-element formed 
with four smoothing triangles in a cross-diagonal 
pattern has been reduced to a four-node quadrilateral by 
exactly enforcing a set of natural edge-wise constraint 
equations. The method, smoothing element analysis, 
provides (^-continuous recovered stress distributions 
based on the minimization of a penalized-discrete-least- 
squares error functional. The error functional involves a 
discrete least- squares term in which discrete finite 
element stresses are compared with continuous recovered 
stresses, a penalty constraint term that enforces C 1 - 
continuity of the recovered stresses, and a curvature- 
control term that ensures stability and robustness of the 
method. The four-node element is shown to provide 
equivalent results with the original triangular element 
while significantly decreasing computational cost. The 
recovered (^-continuous stress distributions have been 
employed in a posteriori error estimation in NASA’s 
COMET-AR finite element code, thus enabling efficient 
adaptive mesh refinement solutions, including those for 
built-up aerospace structures. 
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